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Abstract — Greedy algorithm are in widespread use for sparse 
recovery because of its efficiency. But some evident flaws exists 
in most popular greedy algoritiims, such as CoSaMP, which 
includes unreasonable demands on prior knowledge of target 
signal and excessive sensitivity to random noise. A new greedy 
algorithm called AMOP is proposed in this paper to overcome 
these obstacles. Unlike CoSaMP, AMOP can extract necessary 
information of target signal from sample data adaptively and 
operate normally with little prior knowledge. The recovery error 
of AMOP is well controlled when random noise is presented 
and fades away along with increase of SNR. Moreover, AMOP 
has good robustness on detailed setting of target signal and less 
dependence on structure of measurement matrix. The validity of 
AMOP is verified by theoretical derivation. Extensive simulation 
experiment is performed to illustrate the advantages of AMOP 
over CoSaMP in many respects. AMOP is a good candidate of 
practical greedy algorithm in various applications of Compressed 
Sensing. 

Index Terms — Compressed Sensing, Greedy algorithm. Sparse 
Recovery, Noisy Environment. 



I. Introduction 

S Parse signal recovery problem is the reconstruction of 
such signals with characteristic of "Sparsity" from a set 
of nonadaptive linear measurements. It has great potential of 
application on various engineering fields such as coding and 
information theory, signal processing, machine learning and 
others (See papers in website [1] and reference herein). Sparse 
signals contains much less information than their ambient 
dimension suggests. Most of entries in its vector representation 
are zero (or negligible). So it is possible to reconstruct original 
signal approximately or even accurately using only a small 
number of linear measurements. The measurements are of the 
form where $ is some mxN measurement matrix, m <C 
N and x is original signal. Clearly the process of sparse signal 
recovery could be formulated as solving undetermined linear 
algebraic equation y = $a;, where y is measurement data. It is 
well-known that undetermined equation has infinite solutions 
in general. But if we focus on "Sparse" solution only, situation 
will be different. Although the operation for finding the most 
sparse solution of undetermined linear equation is NP-Hard 
commonly[2], theoretical work in compressed sensing has 
shown that for certain kinds of measurement matrices, it is 
possible when the number of measurements m is nearly linear 
in the sparsity of original signal [3] [4]. 

The two major algorithmic approaches to sparse signal 
recovery are based on Li -minimization and on greedy meth- 
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ods (Matching Pursuit). Finding the most sparse solution of 
undetermined linear equation is a Lq optimization problem: 



.||a;||o, s.t \\<S>x - y\\2 < S, 



(1) 



It could be solved by Li relaxation for some measurement 
matrices $[5]. That is, solving ([TJ is equivalent to solving the 
following Li optimization problem 



1, s.t y||2 < (5, 



(2) 



where ||x||i — \xi\ + ••• + \xn\ for x G C". Recently, 
more stronger sufficient condition called Restricted Isometry 
Property (RIP) on measurement matrix $ to guarantee the 
equivalence of ([T]i and (|2]i was also proposed [6]. It was widely 
accepted that Li-minimization ^ was normal path to com- 
plete sparse signal recovery. (|2]i is essentially a linear program- 
ming problem and technique of convex optimization could 
be utilized to solve it effectively [7]. The Li -minimization 
method provides uniform guarantees for sparse recovery. Once 
the measurement matrix <& satisfies the restricted isometry 
condition, this method works correctly for all sparse signals 
X. The method is also stable, so it works for non-sparse 
signals such as those which are compressible, as well as noisy 
signals. However, the method is based on linear programming, 
and there is no strongly polynomial time algorithm in linear 
programming [8]. But its efficiency was questionable and most 
popular software package of convex programming, such as 
cvx[9], is hard to be used in practical application for its 
low rate of convergence, especially when dimension of target 
signal is large. 

On the other hand. Greedy algorithms are quite fast, both 
theoretically and experimentally. It runs by iterating in general. 
Typically, On each iteration, the inner products of residue 
vector r with the columns of measurement matrix <!> is 
computed and a least squares problem is solved to obtain the 
estimation of original signal on this iteration. It is hoped that 
the convergence of iterating could be ensured and the estimator 
could tend to original signal in fewer steps [10]. 

The typical case of greedy algorithms for sparse recovery in- 
cludes Orthogonal Matching Pursuit (OMP) [11], Regularized 
Orthogonal Matching Pursuit (ROMP) [8] and compressive 
sampling matching pursuit (CoSaMP) [12]. It was shown that 
OMP recovered the sparse signal with high probability and 
had great speed, but it would fail for some sparse signals and 
matrices [13]. The development of ROMP provides a greedy 
algorithm with uniform guarantees for sparse recovery, just 
as that provided by Li -minimization method. Furthermore, 
CoSaMP improves upon these results and provides rigorous 
runtime guarantees. However, there are one disadvantage for 
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these two algorithms. Firstly, sparsity level must be presented 
as prior parameter for algorithms. But it is unknown in most 
practical scenario and must be guessed in advance. Once the 
estimated sparsity level has large difference with actual one, 
the error of algorithm will increase evidently (Although this 
error could be analyzed theoretically [12]). The problem be- 
comes more severe in noisy environment. ROMP and CoSaMP 
can't adapt their rurming process to noise condition when 
noise is presented. Actually, noise is inevitable in engineering 
problem and target signal of our recovery algorithms is always 
buried in it. There exist few "Pure sparse" signal in real world. 
Because the dimension of target signal is unknown, it is always 
given with some margin to avoid possible missing. Hence it is 
hard to extract target signal without including certain amount 
of noise. This not only has influence on accuracy of algorithm, 
but also reduce the speed of convergence for algorithms. In 
fact, some calculation is carried through to estimate noise, 
however, which is useless at all. 

A new greedy algorithm for sparse recovery is presented 
in this paper. Compared with ROMP and CoSaMP, our new 
algorithm need no any prior information on sparsity level of 
target signal. Furthermore, it is a kind of "adaptive" algorithm 
which can inspect the existence of noise and adjust the halt- 
ing condition automatically based on detailed state of noise. 
Besides that, it has uniformly guarantee and good efficiency, 
just as ROMP or CoSaMP. Hence our new algorithm is a 
better choice when signal with unknown sparsity level is 
to be extracted (such as compressible signal) under noisy 
background. This paper is organized as follows: In Section 
2 we introduce our new algorithm. Section 3 describes some 
consequences of the restricted isometry property that pervade 
our analysis. The convergence of theorem is also established 
for sparse signals in Sections 3. Practical consideration for 
algorithm implementation is provided in Section 4. Empirical 
performance and some numerical experiment is described in 
Section 5. Finally, Section 7 presents overall conclusion. 

II. Description of new algorithm 
A. Motivation 

The most difficult and important part of signal reconstruc- 
tion is to identify the locations of the components in the 
target signal. The common approach adopted by most greedy 
algorithms is "Local Approximating", that is, computing the 
inner products between measurement vector y and columns 
of measurement matrix We will obtain observation vector 
u = and use u as "Local Approximation" (or "Proxy") 
of target signal x. Note that $ is a dictionary and v is sparse, 
so y has a sparse representation with respect to the dictionary. 
It is reasonable that only a few entries of u are remarkable, 
which imply the locations of the components of x, and most of 
its entries are comparatively small. Of course, the precondition 
for argument above is that $ must satisfy some condition such 
as RIP. Intuitively, given sparsity level n of x, every n columns 
form approximately an orthonormal system. Therefore, every 
n coordinates of the observation vector u look like correlations 
of the measurement vector y with the orthonormal basis and 
therefore are close in the EucUdean norm to the corresponding 
coefficients of x. 



Popular greedy algorithms, including OMP, ROMP and 
CoSaMP, pay much attention to observation vector u and build 
their estimator of location of components in x based on u. 
OMP uses one biggest coordinate of u. It is argued that using 
only the biggest couldn't provides uniformly guarantee. So 
ROMP makes use of the n biggest coordinates of u, rather 
than just biggest one, and take a further step of regularization 
to improve the performance of algorithm. It should be noted 
that sparsity level n is always unknown. CoSaMP employs 
more coordinates of u, the 2n biggest, to avoid the possible 
leakage of component in x. But n must be guessed to be input 
in ROMP or CoSaMP as important parameter. If guessed n is 
smaller than its true value, correct result can't be found; On 
the other hand, if guessed n is set to a very large value (maybe 
much larger than true value) to ensure that all of entries of x 
will enter the view of algorithms, certain amount of noise will 
presented in our calculation inevitably. How can we make a 
good guess for sparsity level n without any knowledge on its 
true value? 

B. AMOP Algorithm 

Our new approach, named Adaptive Orthogonal Matching 
Pursuit (AMOP), chooses appropriate number of biggest en- 
tries of observation vector u by studying the fine feature of 
u. At each step, u is computed by w = where r is the 
residue vector of last step. Unlike other algorithms, AMOP 
determines the estimation for n by analyzing the trend of 
entries in u arranged by descend order of their amplitude. That 
is, relative amplitude difference of adjacent elements in above 
queue is calculated and a threshold is set. Once the relative 
ampUtude difference between fcth and (fc -|- l)th element in 
ordered queue of entries in u exceed threshold, k will be 
chosen as estimation of n. 

Detailed description of AMOP is proposed as follows: 

Algorithm 1 (Adaptive Orthogonal Matching Pursuit): 
Input: Measurement matrix $ e R'^^^ , Measurement 
vector y e R^, Threshold T, e and K. 

1 Let r = t/, n = 0. 

2 Calculate u = ^"r and |u| = (|wi|, • • • , |uAr|). 

3 Arrange |m| by descend order to obtain 



l"U= (l'"l[l])|w|[2],---,|w|[jv]), 

4 Determine index k as follows, let /? = 1, 

kl[fc+i] - l"l[fc] 



(3) 



k = min li £ {1 



,N} 



m 



>T*(3 



5 ifk>Kand(}< 0.1, set k = K; else (3 = /3* 0.9, goto 
step (4); 

6 Update the set of indices by f2 = U {[1], • • • , [A;]}. 

7 Solve least square problem 

mm\\^\ax-y\\2, 

X 

where is a submatrix of $ composed of its columns 
with index in f2. 

8 Calculate the residue vector of r = y — ^x. 

9 If |r"|/|y| < e, output x and O, stop; else go to step (2). 
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For the descend order of queue (O, we have 



As input, the AOMP algorithm requires two adjustable 
parameter T and e besides matrix $ and measurement vector y. 
But it doesn't need sparsity level of target signal x anymore, 
unlike ROMP and CoSaMP. It can be extracted incidentally 
along with running of algorithm. Furthermore, the number of 
components selected in the step (4) is determined by algorithm 
itself automatically. It is easy to understand that this number 
is critical for performance of algorithm. Any manual setting 
wiU introduce extra error when mismatch between prescribed 
value and actual situation of data exists. So it is very necessary 
to let greedy algorithm of sparse recovery be adaptive, just as 
AOMP 

Step (5) should be noted that it give AMOP algorithm 
more flexibility and stability. If threshold T is set too large 
so that too much coordinates was selected in one iteration, 
algorithm is prone to degrade or crash. For this we build a 
upper bound in step (5) to prevent the crazy growing of number 
of chosen components. If this bound is exceed, threshold T 
will be adjusted to smaller value to increase the possibility of 
components in \u\ in step (3) to satisfy the condition in step 
(4). The importance of step (5) is also illustrated in following 
section on analysis of convergence. 

III. Theoretical Analysis of Algorithmic 
Performance 

There are two kinds of iterative invariant of greedy algo- 
rithm for sparse recovery deduced in the convergence analysis 
for ROMP and CoSaMP respectively. As to CoSaMP, assume 
the sparsity level s is preliminary and 

V ^ \\x - Xsh + -^\\x ~ Xs\\l + ||e|l2, 

where Xg is s-sparse approximation for x and e is additiv 
noise, the following assertion could be proved [12]. 



< lla; 



(4 



where is result of pruning step in CoSaMP. Hence it i 
forced to be s-sparse. So this kind of iterative invariant is nc 
suitable for analysis of AMOP because the sparsity level o 
intermediate result at each step in AMOP isn't fixed. Howevei 
the iterative invariant in ROMP simply concerns with th 
percentage of the coordinates in the newly selected set belon: 
to the support of target signal x. It is argued that the ratii 
above isn't lower than 50% with the help of regularizatio 
step [8]. We will follow the idea in [8] to derive our result oi 
convergence of AMOP. 

A. Localization of Energy 

By induction on the iteration of AOMP, we study the gain in 
each iteration. Losing no generality, suppose sparsity level of 
target signal x be S, and k coordinates is selected eventually 
in this iteration. Then its percentage of energy of first k 
components in queue (O is 



P = 



Sn=l \y\'[n\ 



X]n=l l^l^n] + Sn=fc+1 l^lfn 



(5) 



P > 



> 



[k] 



T.l=MU + iK-k)\y\ 



[k] 



[fe] 



E:=oil-T)-^^'' + iK-k){l-T)^ 
k 



(6) 



VX^^^: + {K-k){i-Tr 

It is easily seen that |6] achieves its minimum at fc = 1, that is 



p — 

^ vn.%n. — 



1 



I + [K -Tf 



(7) 



Here fc = 1 means only the largest coordinates was chosen, 
which is just the choice of OMP. So OMP is a special case 
of AMOP. According to [8] Lemma 3.6, a large portion of 
energy of unidentified part of target signal would be locked 
by queue (O and certain amount of energy would be reserved 
by "regularization step" in ROMP by [8] lemma 3.8. In AMOP, 
the "regularization step" is replaced by choosing k largest 
coordinates. So more energy would be identified in AMOP 
than in OMP because more than one coordinates would be 
chosen in AMOP generally. The ability of locking uncovered 
energy of target signal for AMOP and ROMP is compared in 
FiglD 




10 
Sparsity 



Fig. I . Ability of Locking Uncovered Energy 

It is shown that ability of AMOP is superior to that of 
ROMP when sparsity of target signal is small. Although the 
advantage becomes vague when sparsity grows, AMOP still 
has relatively good performance considering more coordinates 
would be chosen and percentage of actual identified energy 
would be larger than 
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B. Getting the "Correct" Support 

In [14] and [15], the correctness of support of solution 
for OGA algorithm under different noise scenario were ana- 
lyzed. Mutual coherence of matrix $, overcomplete dictionary 
system, was the key parameter for performance of greedy 
algorithm in discussion therein. Here we will give analogous 
results for AMOP under noisy condition based on RIC (Re- 
stricted Isometry Constant)[3] of $. 

In practice the noise is unavoidable and it is always assumed 
that ideal noiseless signal y'^ has sparse representation = 
^x'^, where the support of is very small. What we can 
observe is noisy version y — y'^ + n where ||n||2 < e. Suppose 
a;° be solution of 



mm||a:||o, 



s.t <^x°^y°, 



be final result of AMOP, and 



S'^ = supp(a::"), S — supp(x'*) 



(8) 



(9) 



We argue that 5C5° under certain conditions on value dis- 
tribution of target signal. That is, the correctness of support 
of solution for AMOP can be guaranteed even in noisy 
environment. 

Proposition 1: if target signal {a;^} satisfy 



> 2(e- 



— max kr,- 1 
2 ~^co " ' 



(10) 



jes" " ■'" " 1 - (5k V 2 jes° 

here matrix $ is supposed to has X-order restricted isometry 
constant 5k, K is twice of sparsity level of target signal x^ . 
Then iSC5° holds throughout the iteration process of AMOP 
unless all the coordinates in were chosen. 
Proof: 

We proceed by induction. SQS^ holds at beginning of step 
1 in AMOP initially for 5 = 0. Assume it is true at beginning 
of step 1 in given iteration, we prove it is still true at beginning 
of step 1 in next iteration. Consider case of S<zS^ , we have 

r = 'i'sxs-y, (11) 

It is trivial that 

$52^5 + $SO\s4o\5 - y" = *5o4" - 2/" = 0, (12) 

because xs is the solution of least square optimization on step 
7 in AMOP, 



Multiply ($^$5) on two side of[ 



(13) 



0, 

(14) 



we have 



= ^SXS -y- ^SXs + 



s"\s 



S°\SXso\s 



(15) 



where / is the identity matrix. 

For j ^ S, the norm of (j)Jr is estimated and bounded. 
Firstly, because — / is the projection matrix 



of orthogonal complement of subspace spanned by $5, its 
2-norm is 1. So 



< 



^(<i>5($^<i>5)-'<fl-/)(y~2/°)||2 
M2ms{'i>^'i>sr''^l-m2\\{y-y°)\\2 



1 * 1 * e 



(16) 



Secondly, 



*5 («'s $5 ) " $5 XS^S" \sh 
>5||2||($5*5)"'||2||$5*50\5l|2|l4o\5ll2(17) 

Because tt{5U5°} < 25* and K = 25*, According to [12], 
proposition 3.2, we have 



and 



l*S*50\<sl|2 < Sk 



'$s||2 < Sk 



(18) 



(19) 



On the other hand, according to definition of RIC, we obtain 



^ <ll(*?<i'5)-^l|2< 



OK 



'K 



(20) 



Hence 



1 1 J $5 ) - 1 $50 \sx% \5 1 1 2 



< 



SI 



Thirdly, by analogous deduction, if j ^ 5°, 

Il0j$5«\54«\5ll2 < SK\\x%o\sh 

Summarize the results above, we have for j ^ 5°, 



(21) 



(22) 



\-S 



K 



5k)\\x% 



S0\Sl|2, 



< 



< 



3a- 



1-5 



K 



ll4o\5l|2, 



< e + 



^ max II4II2 
— max ||a;°ll2, 



(23) 



Lower bound for ||0"fr||2 is considered similarly. For j £ 



> 
> 
> 



+ 0j$5«\(5u{i})4o\(5ub})ll2, 

X^jh - \\4'J'^S°\iSU{j})X%a\^su{j})\ 
X^h - '^i<'ll4«\(SUb})ll2 

a;°||2 - SK\\x%,^sh 



(24) 



Hence 

4.T 



Ujrh > \\x%-6k\\x%\s\\2 

5 k „ 



^ ll4o\5ll2 



1~5k 



= \\x% 
> \\x% 



X 



.0 



1 - 5 k 
5 k IK 



5"\Sll2 



1-5k 



2max||.,.||2, 



(25) 



5 



if target signal {5°} satisfy ( fTOl i. we have 



min \\(h, r U > max (A,- 



^2, 



(26) 



That is to say, 5 C 5° will hold throughout the iteration 
process of AMOP unless all the coordinates in were 
chosen. ■ 
It is argued that value distribution of target signal, power 
of noise and RIC of matrix $ are all critical to performance 
of greedy algorithm. The proposition above gives a general 
condition for correctness of support of solutions for a large 
class of greedy algorithms (Not just AMOP) which use inner 
product between residue r and dictionary vectors (columns 
of $) to obtain information of support of target signal. It 
seems that condition ( fTOl i is too restricted. But it is easy to 
see from ( fTOl i that "Dynamic Scope" of target signal (that 
is, the norm difference between the elements with maximal 
and minimal norm) depends on RIC 5k of matrix $ and 
noise power e. Consider the requirement on RIC in ROMP, 
which is 0.03/-\/log(s) with s is sparsity level of target signal 
according to [8], Theorem 3.1, we write ( fTOl l as 



0.06^/i 



0.03 



2e, 



(27) 



with maximum is normalized to 1. It is depicted in Fig|2] 
for noise level is 0.1(SNR is 20dB). When sparsity level is 
small, target signals with considerable 'Dynamic Scope' are 
guaranteed to have good performance in greedy algorithms. 
The restriction on 'Dynamic Scope' of target signal becomes 
tighter gently when sparsity level increases. The actual number 
of chosen coordinates in iteration step of AMOP in practical 
scenario is smaller than sparsity level in general. So AMOP 
could choose correct coordinates in most cases. This assertion 
will be illustrated further in numerical experiments. 
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IV. Practical Consideration For Algorithm 
Implementation 

A. Least Square via Orthogonalization 

For efficient implementation of AMOP, it is important to 
design a appropriate computational scheme with low com- 
plexity and good numerical stability to calculate the solution 
of least square problem in step 7 in AMOP. It should be 
noticed that AMOP is incremental, that is, in each iteration 
some new coordinates were selected and none of previously 
chosen coordinates was excluded, unlike CoSaMP. So it is 
possible to construct a recursive algorithm to solve the least 
square problem. 

Assume on the first iteration several coordinates were cho- 
sen and denote the set of corresponding columns of matrix $ 
as Ai. Then observation vector y would be linearly approx- 
imating with vectors in Ai and coefficients xiare computed 
by solving the following equation 



A^A^xi = A^y, 



(28) 



ordinary solver with good numerical performance such as QR 
decomposition or Singular value decomposition could be used 
to compute xi. From geometrical point of view, calculation 
of xi is equivalent to project y onto subspace spanned by Ai. 
We wrote it intuitively as 



xi = y\Ai, 



(29) 



On the next iteration, a set A2 of columns of matrix $ with 
respect to newly chosen coordinates would be added to least 
square regression. The projection became 



X2 = v\{Ai,A2}, 



(30) 



It is well-known that orthogonalization could simplify the 
calculation for projection onto subspace spanned by two 
mutually orthogonal subspace could be regarded as sum of 
projections on each one. So A2 was written as 



A2 = A\®A^ 



(31) 



where A\ was projection of A2 onto Ai and A2 is orthogonal 
to ^1. Hence 



X2 = y\{A^,A2} = y\{A,.Al} 
= y\Ai+y\Al^xi+y\Al, 



(32) 



This could be accomplished with Gram-Schmidt orthogonal- 
ization procedure. Without loss of generality, suppose 



A2 



then 



A1UA2 = {01 >2, 

using following procedure 



(33) 



Uk 



k-l 

E 



( ^tn 1 Ufa ) 



(34) 
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where (•) denotes the inner product of vectors in Euclidean 
space. We can obtain 

Bi = iUuU2,---,Uk), 

B2 = {Uk+l,Uk+2,---,Un). 

The projection in (|32] | could be written as 



y\Al = y\B2 = y\Uk+i 



-y\Un 



(35) 



The numerical stability of Gram-Schmidt orthogonalization 
procedure could be improved further [16]. Instead of com- 
puting the vector C/fc as (|34] i. it is computed as 

Vk - 777— TTT^l' 



(1) 



{Ui,Ui 

(1) {ui'\u2) 

{U2,U2) 
Ak-3) 



U2 



= u, 



(/c-2). 



Uk 



(fe-2) 



{ur'^^Uk 



_^(k-2) 



(fe- 



1) Tjik- 



-u. 



(fc-1) 



(36) 



This approach gives the same result as the original formula 
in exact arithmetic, but it introduces smaller errors in finite- 
precision arithmetic. 

There are one points worth mentioning. Although other 
orthogonalization algorithms such as Householder transforma- 
tions or Givens rotations are more stable than the stabilized 
Gram-Schmidt process, they produce all the orthogonalized 
vectors only at the end. On the contrary, the Gram-Schmidt 
procedure produces the jth orthogonalized vector after the 
jth iteration and this makes it the only choice for iterative 
algorithm like AMOR 

B. Resource Requirements 

AMOP was designed to be a practical method for sparse 
signal recovery. The main barrier for algorithm efficiency is 
least square problem in step 7 of AMOP. Using recursive or- 
thogonalization procedure above could mitigate computational 
burden of AMOP dramatically. Furthermore, the orthogonal- 
ization technique has the additional advantage that they only 
interact with the matrix $ through its action on vectors. In 
fact, it only concern with the inner products and additions of 
columns of matrix $. It follows that the algorithm performs 
better when this sampling matrix has a fast matrix-vector 
multiply, such as on parallel computational platforms. On the 
other hand, less memory consumption is another advantage of 
recursive orthogonalization based least square. In fact, direct 
method such as SVD and QR have storage cost O(fcm), where 
k is the number of chosen coordinates in each iteration and 
m is row number of matrix $. It is too huge for large scale 
problems. But for AMOP, only one vector need be put in 
memory in recursive orthogonalization calculation and storage 
cost is 0{m). It is more suitable for implemented with VLSI 
circuit. 

We estimate the time complexity of main steps in AMOP 
as follows: 



TABLE I 
Time Complexity of AMOP 



Step 


Standard 


Fast 


2 


0(mN) 


0{C) 


3 


0{N log N) 


0{N log N) 


4&5 


0{K) 


0{K) 


7 


0{mK) 


OiK) 


Total 


0(niN) 


0(C) 



Step 2: In this step, the inner products of residue vector 
and columns of matrix $ is computed as proxy for sup- 
port of target signal. Its cost is bounded by matrix-vector 
multiply ^^r, which is 0{mN) with standard multiply 
operation or 0{C) for fast matrix- vector multiply. 
Step 3: According to standard textbook on algorithms 
[17], the expected time for selecting largest s entries 
in vector with dimension N is 0{KN). Using efficient 
schemes such as Quicksort or HeapSort, a fully sorting 
of vector could be completed with expected time cost 
0{N \ogN) and largest s entries could be selected di- 
rectly, which is faster n some situation. 
Step 4 & 5: Certain amount of support of target signal 
would be identified in these two steps. Although some- 
times the threshold needs to be adjusted according to step 
5 and several cycles of operations may be necessary, the 
total cost is still 0{K). 

Step 7: The main advantage of AMOP is recursive 
orthogonalization based implementation of least square 
problem. Inner products of vectors are involved in orthog- 
onalization and occupy much of computational resource 
which can be implemented efficiently by matrix-vector 
multiply. The cost is 0{mK) with standard multiply 
operation or 0{£) for fast matrix- vector multiply. 

Table 1 summarizes the discussion above in standard mul- 
tiply operation and fast matrix-vector multiply with cost L, 
respectively. 

Storage cost for AMOP is also considered for showing its 
practicability. Aside from the storage required by the sampling 
matrix, AMOP algorithm constructs only one vector of length 
N as the signal proxy. The sample vectors have length m, so 
they require 0{m) storage. The signal approximations require 
at most 0{s) storage. Similarly, the index sets that appear 
require only 0{s) storage. The total storage is at worst 0{N). 

V. Numerical Experiment 

In this section some numerical experiments were conducted 
to illustrate the performance of signal recovery of AOMP. 
There are three factors to be considered in numerical testing of 
AMOP: construction of $ matrix, value distribution of target 
signal and SNR condition which will be examined in our 
experiments. 
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A. Construction of Matrix $ 

The property of measurement matrix $ is critical to per 
formance of any greedy algorithm for sparse recovery. A 
indicated in section on theoretical analysis. Its RIC has direc 
influence on probability of recovery of algorithms. Here 
several kinds of matrix $ were built and utilized to tes 
the performance of AMOP, including well-known Gaussiai 
random matrix, BernoulU random matrix and random Fourie 
matrix. 

The target signal was set to be flat and no noise wa 
added in, then 500 independent trials were performed. Figur 
|3]|5] depicts the percentage (from the 500 trials) of spars 
flat signals that were reconstructed exactly when m x I" 
Gaussian random matrix was chosen as measurement matri: 
$ . This plot was generated with N = 256 for various level 
of sparsity S. The horizontal axis represents the number o 
measurements m, and the vertical axis represents the exac 
recovery percentage. 

As comparison, recovery percentage of algorithm CoSaMP 
is also given under the same setting. Standard CoSaMP needs 
sparsity of target signal as its important prior knowledge am 
it is widely regarded as one of main drawbacks of CoSaMP. Ii 
our experiment, sparsity of target signal was given to CoSaM) 
as input parameter to guarantee the power of CoSaMP to b 
exploited fully. 
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Fig. 3. Recovery Percentage of Signal with Sparsity 4 with Gaussian Matrix 

It should be noted that even the sparsity of target signal is 
known beforehand (which is impossible in practice), recovery 
percentage of CoSaMP is lower than that of AMOP. Especially 
when sparsity was relatively large, performance of CoSaMP 
degenerated very rapidly. On the contrary, the behavior of 
AMOP was very stable. According to well-known theoretical 
result of Compressed Sensing, for Gaussian random measure- 
ments matrix if row number m, column number TV and 
sparsity S satisfies 
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Fig. 4. Recovery Percentage of Signal with Sparsity 20 with Gaussian Matrix 
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m>CS\og[N), 



(37) 



Fig. 5. Recovery Percentage of Signal with Sparsity 36 with Gaussian Matrix 



where C is a constant independent of 5, then the probability 
of recovery failure is exponentially small [3]. Our experiment 
result indicates that for AMOP, the value of constant C is 
about 2 for Gaussian measurement matrix. 

Figure |6]|8] depicts corresponding result for Bernoulli ran- 
dom measurement matrix, which is analogous to Gaussian 
case. It had been proved that condition ( [37] is also sufficient 
for overwhelming probability of successful recovery for binary 
Bernoulli measurement matrix [18]. It is observed that the 
constant C for AMOP in Bernoulli case is probably the same 
as that in Gaussian case. 

Figure |9][TT] depicts corresponding result for Fourier random 
measurement matrix. Somewhat surprisingly, it is similar to 
that of Gaussian and Bernoulli case. To our knowledge, the 
best known bounds on size of measurements in Fourier case 
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Fig. 6. Recovery Percentage of Signal with Sparsity 4 witli Bernoulli Matrix Fig. 8. Recovery Percentage of Signal witli Sparsity 36 witli Bernoulli Matrix 
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Fig. 7. Recovery Percentage of Signal with Sparsity 20 with Bernoulli Matrix pig. 9. Recovery Percentage of Signal with Sparsity 4 with Fourier Matrix 



is given by [19] 



m>CS{log{N))^ 



(38) 



which is conjectured to be the same as |37] [3] but there exists 
no strict theoretical proof until now. Our experiment result 
verified this conjecture in some extent indirectly. 

B. Value Distribution of Target Signal 

Pure flat signal is rarely seen in practical engineering 
application. So it is necessary to investigate the performance 
of sparse recovery algorithms on non-flat target signal. There 
are two cases to be studied. One is piecewise flat signal 
which is common in various fields of imaging, such as optical, 
microwave and magnetic resonance. The result is depicted in 
Figure [12] to [14] Here the measurement matrix is fixed to 
Gaussian random matrix. 



Recovery percentage of CoSaMP in our experiment isn't 
ideal. Especially when the number of non-zero elements of 
target signal is large, the successful recovery probability of 
CoSaMP becomes more and more ignorable. The capability 
of CoSaMP is doubtful in this setting. On the contrary, the 
behavior of AMOP is very robust when value distribution of 
target signal is changed. Furthermore, the constant C in [37] 
is approximately equal to that in flat signal setting. Although 
it is predicted theoretically that constant C only depends on 
the construction of measurement matrix $, not rely on nature 
of target signal, it is common in practice that the detailed 
value distribution of target signal certainly has influence on 
performance of recovery algorithms. Our experimental result 
indicates that the actual behavior of AMOP coincides with 
theoretical conclusion perfectly. This argument is confirmed by 
result for the other case of non-flat target signal. Here target 
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Fig. 10. Recovery Percentage of Signal with Sparsity 20 with Fourier Matrix Fig. 12. Recovery Percentage of Piecewise Flat Signal with Sparsity 4 
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Fig. 11. Recovery Percentage of Signal with Sparsity 36 with Fourier Matrix Fig. 13. Recovery Percentage of Piecewise Flat Signal with Sparsity 20 



signal is chosen as compressible signal which is frequently 
presented in orthogonal representation of signal, such as 
Discrete Cosine Transform and wavelet transform, and data 
compression. Two kinds of compressible signal are analyzed 
in our experiment, one is exponential signal, 

x{n) ^ C -^a^ , 0<a<l 

the other is polynomial signal 

= C*7ii/P, 0<p<l 

For brevity, only result for exponential signal is depicted in 
[TS] and corresponding curve for CoSaMP is omitted. We also 
performed the same experiment for polynomial compressible 
signals and found the results very similar to those in Figure 

m 



C. Target Signal With Noise 

Random noise was added in target signal to test the per- 
formance of sparse recovery algorithms in noisy environment. 
Measurement matrix is fixed to Fourier random matrix and 
target signal is set to flat. The sparsity of target signal is fixed 
to 20 and size of measurements m is fixed to 200. The relative 
error of AMOP and CoSaMP in Gaussian white noise with 
various level is depicted in Figure [16] 

The capacity of AMOP in noisy environment is satisfactory. 
Relative recovery error could be controlled within 10% when 
SNR is about lOdB. Even when SNR is as low as 5dB, relative 
error of AMOP still could be governed within 20%. Though 
the error curve rise very acutely in very low SNR region, it is 
shown that AMOP works normally in most noisy environment. 
On the other hand, the relative error of CoSaMP keeps on a 
high level when noise is presented. When SNR was increased. 
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Fig. 15. Recovery Percentage of Polynomial Flat Signal 



it didn't exhibit the obvious trend of decreasing. With the 
language of statistics, CoSaMP isn't consistent estimator 

Figure [17] illustrates the advantage of AMOP in noisy 
environment from other viewpoint. Here SNR is fixed to lOdB, 
the error with various size of measurements is plotted. It 
is observed that if number of measurements n is small, the 
performance of AMOP is heavily abnormal. In fact, relative 
error of AMOP is even higher than CoSaMP when n is lower 
than 100. But it falls abruptly when m increases while that 
of CoSaMP remains in narrow range. When m is larger than 
100, the relative error of AMOP tends to be stable. It is well 
controlled with in 10%, which is similar to Figure [16] 



D. Measurement Matrix in STAP 

We would build a measurement matrix $ with spatial- 
temporal basis vectors, which is key component in theory 
and computation of STAP (Space-Time Adaptive Processing), 
to investigate the potential of sparse recovery algorithms 
to be applied in field of modern radar and communication 
engineering. Here $ is set to 

$ = [0,_t(l, 1), • • • , (Ps-t{l, n),---, (Ps-t{m, n)], (39) 

and 

(/Js-tifsJd) = [l,exp(j27r/d),---,exp(j27r(L-l)/d), 
• • • , exp(j27r((7V - 1)/,, + [L - l)fd))f, 

where and denotes Doppler and spatial frequency, 
respectively. Unlike Fourier matrix, the construction of spatial- 
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temporal matrix <i> is more complex. The phase of elements 
in $ composed of two parts, one is contributed by Doppler 
frequency and the other by spatial frequency. It lead to 
following consequence: Different fd and fs could be combined 
to form the same (or approximately equal) phase. In other 
words, strong correlation exists in different column vectors 
in $, which corresponding to various points far away with 
each other on spatial-temporal plane. So Restricted Isometric 
Constant of $ is conjectured to be relatively large. It is 
challenge for sparse recovery algorithm to be feasible whei 
measurement matrix $ is chosen as spatial-temporal matrix. 

Generally speaking, the support ("Position" of non-zer 
elements) of target signal is much important than its detailei 
value. It is indeed true in practical engineering discipline. Fo 
example, in radar STAP processing, sparse recovery is utilizei 
to estimate the energy distribution of clutter and interferenc 
on spatial-temporal plane from sample data directly. Becaus 
echo of clutter and interference is much stronger than rada 
target, the detailed amplitude and phase of clutter and interfer 
ence isn't crucial. As long as the accurate support ("Position" 
of clutter and interference on spatial-temporal plane is foum 
out, we can design the efficient filter to suppress clutte 
and interference effectively. Hence the most important featur 
of sparse recovery algorithm applied to STAP processing i 
detecting the support of target signal with great precision. 

The experiment was performed to test the ability of AMOl 
and CoSaMP to detect signal support. Measurement matrix 
$ was set to 224x900 spatial-temporal matrix as [39] Target 
signal is noise-free flat signal and its sparsity is fixed to 20. 
Figure [18] depicts the average result from 100 trials. 
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Fig. 18. Size of Detected Target Support With No EiTor 

It is clear that AMOP is much superior to CoSaMP when 
applied to STAP processing because it could detect a majority 
of target support with no error while CoSaMP could only 
find very few. But even so, the accuracy of algorithms for 
support, whether AMOP or CoSaMP, can hardly satisfy the 
requirement of practical STAP processors. Due to ultra-low 



Signal Clutter Ratio (generally lower than -50dB), missing 
four or five frequency points on spatial-temporal plane would 
lead to very high false alarm rate and the performance of radar 
would degrade heavily. So we should detect as much target 
support as possible to minimize false alarm rate. If tiny error 
is allowed, the behavior of sparse recovery algorithms becomes 
better, as depicted in Figure [T9l 
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Fig. 19. Size of Detected Target Support With Error 1 

It is easily seen that if the "Position" detected having 
difference 1 with true "Position" of target signal is allowed 
to be counted, the average size of detected target support 
for AMOP increases by about 2 and that of CoSaMP is still 
negligible. It accounted for that some support of target signal 
missed by AMOP wasn't really missed. That is to say, their 
neighborhood, which corresponding to the points adjacent to 
them on spatial-temporal plane, were discovered instead. This 
is a good news for high performance filter to suppress clutter 
could still be designed with AMOP and carefully chosen 
notch, without losing much resolution and SNR. Figure |20] 
depicted the case of error 2. The behavior of AMOP continued 
to be made better and approach the best. It seems that it make 
little sense to enlarge error tolerance further. 

VI. Conclusion 

In this paper, a novel greedy algorithm for sparse recovery, 
called AMOP, was given and examined. Its performance was 
studied by theoretical analysis of simulation experiment. 

The motivation of this algorithm is two obvious drawbacks 
in popular methods, such as CoSaMP: Need of sparsity of 
target signal as prior knowledge, and weak ability of working 
in noisy environment. With well-designed algorithmic steps, 
AMOP can extract the information on sparsity of target signal 
adaptively and sense the nature of target signal automatically. 
It can recover the detail value of target signal with very 
high precision with little prior knowledge. Its validity is 
illustrated by strict deduction. Fine stability of performance 
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Fig. 20. Size of Detected Target Support With EiTor 2 



under random noise is another advantage of AMOR Fur- 
thermore, its robustness for various setting of target signal, 
flat or compressible, and construction of measurement matrix, 
such as spatial-temporal matrix, were also shown by thorough 
numerical experiment. It is argued that AMOP is a excellent 
greedy algorithm for sparse recovery and has great potential 
of widely utilization on signal processing. 
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